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Quantum chromodynamics with two zero mass flavors is expected to exhibit a phase transition 
with 0(4) critical behavior. Fixing the universality class is important for phenomenology and for 
facilitating the extrapolation of simulation data to physical quark mass values. Other groups have 
reported results from lattice QCD simulations with dynamical staggered quarks at Nt = 4, which 
suggest a departure from the expected critical behavior. We have pushed simulations to the largest 
volumes and smallest quark mass to date. Strong discrepancies in critical exponents and the scaling 
equation of state persist. 



I. INTRODUCTION 



It is generally expected that two-flavor QCD undergoes a high temperature chiral-symmetry-restoring phase tran- 
sition at zero quark mass with 0(4) critical behavior Should the axial anomaly disappear simultaneously with 
the phase transition, the Pisarski Wilczek scenario then suggests a fluctuation-driven first order phase transition. 
Verifying these expectations is important for understanding the phenomenology of the transition and for facilitating 
an extrapolation of simulation data to physical quark masses. Since the staggered fermion scheme breaks the anomaly 
explicitly at nonzero lattice spacing, lattice QCD at fixed Nt with staggered fermions, as a statistical system in its 
own right, is similarly expected to exhibit at least 0(2) universality, with 0(4) or a fluctuation-driven first-order 
phase transition emerging in the continuum limit Nt oo. At TVt = 4 the lattice spacing is coarse enough that, if 
there is a critical point at zero quark mass, 0(2) is the only likely option. 

The standard test of universality compares critical exponents. Comparing the critical scaling function itself gives 
further insight. To test for the expected universality we use the standard correspondence between QCD variables 
and 0(N) spin variables, which identifies quark mass niq/T with magnetic field h, inverse gauge coupling Q/g^ with 
temperature T/Tc(0), chiral condensate (iptp) with magnetization M, and the action (plaquette) with the energy 
density. A critical point is expected to occur at zero quark mass and nonzero coupling 6/g^(0). For studies at fixed 
Nt, therefore, we define |] 

h — aniqNt 

t = - 6/gj ex T/T, ~ 1 (1) 
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Critical scaling theory predicts that for small quark masses we have the Fisher scaling relation 



{^^)h-'^' = facDix) = cyfaic^x) (2) 

where x = th~^^^^ , Cx and Cy are scale constants, ^Qcoix) is the critical scaling function for QCD and fcix) is that 
for the appropriate universality class G. Only the scale constants Cx and Cy are adjustable. Some critical exponents 
arc given in Table ^. Outside the Ginzburg scaling region, by definition, there are appreciable nonleading, nonscaling 
contributions to analytic in t and h. In addition to corrections analytic in t and h, there are correction terms 

with subleading exponents, universal and nonanalytic in t and h. The mean-field scaling function is known exactly. 
For the 0(4) scaling function we use results of a numerical simulation [Q. 

There is a similar scaling relation for the energy density. In QCD the energy density (plaquette) is dominated 
by gluon degrees of freedom, which are indirectly affected by the chiral singularity. So apparently there is a much 
larger analytic contribution. Consequently, we have found the plaquette observable much less useful for testing critical 
scaling. Here we concentrate on the scaling of {'ipip). 

In the next section we discuss an analysis of finite size effects, present a determination of some critical exponents, 
and compare our results with the critical scaling function. In the concluding section we suggest reasons for the 
discrepancies observed. A preliminary version of this study was presented at Lattice '97 H. 



II. ANALYSIS OF CRITICAL BEHAVIOR 



Our data set extends an old sample on lattice sizes x Nt with A^t = 4, L = 8 and quark masses aniq — 0.025 and 
0.0125, which was generated with the standard one-plaquette gauge action plus two-flavor staggered fermion action. 
Our new simulations decrease the quark mass to amq = 0.008 and increase the spatial lattice size L to 24 (aspect 
ratio 6). We also reanalyzed old data at Nt = 6, 8, and 12 ||^. The old data, unfortunately, are limited to aspect 
ratio L/Nt = 2. The extent of our iVj = 4 data sample is given in Table ||. Included in this table are values for global 
observables. For equilibration we typically dropped the first 300 molecular dynamics time units of each run. 

Over the range of nonzero quark masses considered, there appears to be no phase transition - only a crossover, as 
illustrated for aruq = 0.008 in Fig. |l|. Evidently, however, the crossover steepens as the lattice volume is increased. 
The crossover, or "pseudo-critical point" is signaled by a peak in a susceptibility for any lattice size. For example the 
mixed plaquette/chiral condensate susceptibility, corresponding to the slope in Fig. ^ 

" 9(6/g2) ^''> 

is plotted in Fig. || for the 16'^ x 4 lattice. Here, as well as in Fig. |], we use multihistogram reweighting to interpolate 
the data from the simulation points and locate the peak. The error analysis was performed with the jackknife method, 
which enables us to obtain reliable error estimates for both the peak height and location. 

The peak location (crossover coupling Q/gp^) shows little variation in lattice size for L > 8. For example the peak 
location in Xmt shifts from 5.2605(10) to 5.2623(6) as L increases from 12 to 24, a scarcely significant change. It 
also shows little variation among the susceptibilities chosen. For example the peak location varies by ±0.001 over the 
susceptibilities considered. In all cases we take the result from the largest volume and assign an error of 0.002. Close 
to the critical point the peak position occurs at a fixed value of the scaling variable x = Xpc, so we have the scaling 
relation ||] 

6/5pc = 6/5pc(0)2 + Xp.iamqNtY"". (4) 

Shown in Fig. ^ is the trajectory of the pseudocritical point, fitted to both 0(4) and mean-field predictions. Both fits 
are good. An 0(2) fit would do equally well. Such agreement was first found by Karsch and Laermann and inspired 
hope that the simulations had entered the scaling region ||] . 

Problems with scaling were uncovered in studies at larger volumes j7|-|^. In Figs. || and || we plot the peak height 
of the Xmt susceptibility for two fixed quark masses. The increase in peak height with increasing volume reflects 
the steepening trend seen, for example, in Fig. |^. It is necessary to extrapolate to infinite volume at each quark 
mass before checking scaling. We start by assuming the conventional scenario, in which the critical point occurs at 
arUq = 0. Then at nonzero mass, the susceptibility has a finite limit at large volume. We make an ad hoc choice for 
an extrapolation formula with the result shown in Figs. ^ and ^. 

xZT{L) = xZT{^)+blL. (5) 
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The am = 0.0125 data covers the largest range of lattice sizes. Varying the inverse power of L from 1/2 to 1 to 2 
gives a slight preference for 1/L at this mass. Given the uncertainties in the values themselves, we feel it is safe to 
use any of these extrapolations, and we have chosen 1/L for all masses. 

The extrapolated peak height of the Xmt susceptibility is expected to scale with decreasing quark mass. For this 
susceptibility, the expected scaling relation is 

xZT-iam,r-'^/^'. (6) 

We compare this prediction with results from our analysis in Fig. ^ Also shown are similarly extrapolated JLQCD 
values If we include all points in the fit, the scaling exponent is —1.08(8), compared with an 0(4) prediction 
of —0.33 — a clear disagreement, corroborating results of the JLQCD and Bielefeld groups 0-^. However, it is 
evident in Fig. || that at the three lightest masses this observable alone does not exclude 0(4). To test sensitivity 
to our extrapolation formula, we carried out the same analysis, replacing 1/L by l/\fL and 1/L^. The resulting 
scaling exponents are —1.24(11) and —0.94(5), respectively, still clearly at variance with 0(4) over the full mass range 
studied. 

A similar fit of the plaquctte susceptibility, Xtt — d{0) /d{6/g'^), also including 1/L-extrapolated JLQCD results, 
yields a scaling exponent of -0.78(7), while the 0(4) prediction is l/S - 2/((5/3) + 1 = 0.13. The 1/VX and l/L^ 
extrapolations give —0.95(10) and —0.64(5), respectively. 

Because the crossover steepens so much with increasing lattice volume and small quark mass, it is worthwhile looking 
for evidence for two-phase metastability, signaling a first-order phase transition. Figure ^ shows the simulation time 
histories of our runs at arUq = 0.008. While we certainly see long correlation times, we see no evidence for a first 
order transition in these histories. In Fig. || we show time histories from hot and cold starts at arUq = 0.0125 (two 
hot and two cold starts) and am, = 0.008 at values of 6/5^ near the peaks of the susceptibilities. Again, there is no 
evidence for a first order transition. 

We turn next to an analysis of the critical scaling function, given by Eq (||). Here, again, we assume that we are in 
the scaling region. The analysis then depends on which critical exponents we adopt. Using 0(4) critical exponents 
from Kanaya and Kaya |ic| ], we construct /qcd(2;) and compare with the scaling function fo{4){x) for 0(4) Q in 
Fig. ^. Essentially all of the data lie at positive x, which permits a log-log plot. Vertical and horizontal displacements 
of the log- log scaling curves correspond to adjusting and Cy. Clearly, no such displacement would result in good 
agreement. The newer data are plotted with octagons and squares. We observe: (1) The QCD curve falls with 
increasing steepness as the quark mass is decreased. Since the slope of the curve at the crossover gives the peak 
height of the Xmt susceptibility, the disagreement there is consistent with the observed lack of scaling of the peak 
height itself. (2) The new data at larger volume and smaller quark mass show generally worse agreement with the 
0(4) scaling curve. (3) The crossover regions, indicated in the QCD results by line segments and in the 0(4) scaling 
function by a dashed line, are far from agreeing. 

We show a similar comparison of the QCD scaling function with the mean-field prediction in Fig. |l^. Again the 
disagreement is significant. Although we have not measured the 0{2) scaling function, so cannot make a direct 
comparison, given the close similarity of the critical exponents with 0(4), we do not expect any improvement with 
that choice. 

We conclude that if the Nt = A theory falls in the 0(2) or 0(4) universality class, simulations at present masses do 
not reach the critical scaling region. Furthermore, as the quark mass is decreased over the present range, disagreement 
with scaling predictions worsens, offering little hope that we might be getting closer. 

A similar analysis at larger Nt is shown in Fig. |ll|. Perhaps there is improvement with increasing Nt. However, our 
Nt ~ 12 sample includes data only at a single quark mass, making it the weakest test. Furthermore, for Nt > A we 
have no results for L > 2Nt, where we first encountered difficulties at Nt — 4. 

III. DISCUSSION AND SPECULATIONS 

We have seen that new simulations at smaller quark mass and larger volume at A^t = 4 have raised doubts about 
the extent of the previously observed agreement between QCD and 0(4) The conventional staggered fermion 

action with the conventional choice of scaling variables does not show good agreement with the 0(4) or mean field 
scaling functions at present quark masses and temperatures. (Wilson quarks with an improved gauge action seem to 
behave very differently 0.) 

Believers in the conventional sigma model scenario could argue that the critical region is attained only when tt 
and cr correlation lengths are considerably greater than l/Tc- Only in that case is the reduction of QCD to a three- 
dimensional sigma model well justified. Here, typically, these correlation lengths are smaller than 1/Tc. Still, the 
observed worsening of the agreement with decreasing quark mass is disturbing. 
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Recent results from simulations of the conventional Nt = A staggered fermion action, augmented by a four-fermi 
term ("chiral QCD") permit another speculation jlj]. With the additional four-fermi interaction, Kogut, Lagae, and 
Sinclair are able to carry out simulations at precisely zero quark mass. They find evidence for a first order phase 
transition at small four-fermi coupling. A nearby first order phase transition could spoil the approach to the critical 
point. Indeed, one cannot, therefore, rule out the possibility that the first order phase transition extends to zero 
four-fcrmion coupling for a small range of quark masses below the reach of our simulations. In this case, one expects 
a critical end-point at a nonzero quark mass mqc in the Ising or mean- field universality class. At iVt = 6 the same 
group finds evidence for a crossover instead of a first-order phase transition [p^ . Thus, one may speculate that 
the conventional one-plaquette, staggered fermion action at A^t = 4 is plagued by lattice artifacts large enough to 
obliterate the expected aruq = critical point, but these artifacts diminish at higher Nt- 
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TABLE I. Some critical exponents in three dimensions 

yt yh S p 

MF 1.5 2.25 3 0.5 

0(2) 1.495 2.484 4.81 0.3455 

0(4) 1.337 2.487 4.851 0.3836 

Z{2) 1.61 2.5 5.0 0.31 
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TABLE II. Parameters in A^t = 4 data set and two global observables. Run and step lengths are in molecular dynamics time 
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FIG. 1. Order parameter vs inverse gauge coupling for various lattice sizes for am, = 0.008. Curves show results from 
reweighting the data sample. 
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FIG. 2. The mixed (V''!/')~pla'quette susceptibility Xmt as a function of inverse coupling for aniq = 0.008 on a 16^ x 4 
lattice. Curves show results from reweighting together with one-standard-deviation bootstrap errors. 
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FIG. 3. Trajectory of the pseudocritical point G/gin as a function of quark mass in units of temperature mq/T for Nt — 4. 
Crosses indicate points from Karsch and Laermann pi. Also shown are fits to both mean-field (upper curve) and 0(4) (lower 
curve) scaling predictions. 
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FIG. 6. Scaling analysis of extrapolated peak height for the mixed (?/)^)-plaquette susceptibility Xmt for Nt — 4. Results 
from Ref. M are obtained by similar infinite volume extrapolation. Also shown is the 0(4) scaling prediction. 
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FIG. 7. Simulation time history of the order parameter {tf>4') near the crossover for the largest volume 24^ x 4. 
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FIG. 8. Same as Fig. | showing evolution from hot and cold starts at or near the crossover. In each case the initially upper 
(lower) trace follow a cold (hot) start. 
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FIG. 9. Scaling test at A^t = 4 of the order parameter {^'4)), based on 0(4) critical exponents. Shown for comparison is the 
0(4) scaling function from Ref. The crossover region is indicated by line segments in the data and a dashed line in the 
0(4) scaling function. 
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FIG. 10. Same as Fig. but with mean field exponents and the mean field scaling function. 
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FIG. 11. Same as Fig. | but for Nt = 4, 6, 8, and 12. 
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